if !d.name eq 'PS' then begin
    device,xsize=18,ysize=12,yoffset=3
    !p.thick=2 & !x.thick=1 & !y.thick=1
end

@eu_setup
@rhoprof

ns=110D
ng=100D
rsinth1=fltarr(m,l)
d_fr_dphi=fltarr(n,m,l,ns)
d_fth_dphi=fltarr(n,m,l,ns)
dsin_fphi_dth=fltarr(n,m,l,ns)
d_fr_dth=fltarr(n,m,l,ns)
dr_fth_dr=fltarr(n,m,l,ns)
dr_fphi_dr=fltarr(n,m,l,ns)
;
fphi=bx
fth=by
fr=bz
;
; curl components
cfr=fltarr(n,m,l,ns)
cfth=fltarr(n,m,l,ns)
cfphi=fltarr(n,m,l,ns)

;TOTAL HELICITY COMPONENT

for it=0,ns-1 do begin
print,'Computing total helicity',it
 for k=0, l-1 do begin
  for i=0,n-1 do begin
   dsin_fphi_dth[i,*,k,it]=xder_curl(reform(smooth(sin(theta)*fphi[i,*,k,it],3)),theta)       
   d_fr_dth[i,*,k,it]=xder_curl(reform(smooth(fr[i,*,k,it],3)),theta)
  endfor
 endfor

 for j=0,m-1 do begin
  for i=0,n-1 do begin
   dr_fth_dr[i,j,*,it]=xder_curl(smooth(r[*]*fth[i,j,*,it],3),r)   
   dr_fphi_dr[i,j,*,it]=xder_curl(smooth(r[*]*fphi[i,j,*,it],3),r)
  endfor
 endfor
;

 for k=0, l-1 do begin
  for j=1,m-1 do begin
   rsinth1[j,k]=1./(r[k]*sin(theta[j]))
  endfor
 endfor
 rsinth1[0,*]=0.
 rsinth1[m-1,*]=0.


 for k=0, l-1 do begin
  for j=1,m-1 do begin
   for i=0,n-1 do begin
     cfr[i,j,k,it]=rsinth1[j,k]*dsin_fphi_dth[i,j,k,it]
     cfth[i,j,k,it]=rsinth1[j,k]*d_fr_dphi[i,j,k,it]-dr_fphi_dr[i,j,k,it]/r[k]
     cfphi[i,j,k,it]=dr_fth_dr[i,j,k,it]/r[k]-d_fr_dth[i,j,k,it]/r[k]
   endfor
  endfor
 endfor

endfor

mht=total((cfr*fr + cfth*fth + cfphi*fphi),1)/double(n)
mhtm=total(mht[*,*,ns-ng:ns-1],3)/double(ng)

; CURL OF THE MEAN
print,'Computing mean helicity',ns
ffphim=total(fphi,1)/double(n)
ffthm =total(fth,1)/double(n)
ffrm  =total(fr,1)/double(n)

fphim=total(ffphim[*,*,ns-ng:ns-1],3)/double(ng)
fthm=total(ffthm[*,*,ns-ng:ns-1],3)/double(ng)
frm=total(ffrm[*,*,ns-ng:ns-1],3)/double(ng)

d_fr_dphim=fltarr(m,l)
d_fth_dphim=fltarr(m,l)
dsin_fphi_dthm=fltarr(m,l)
d_fr_dthm=fltarr(m,l)
dr_fth_drm=fltarr(m,l)
dr_fphi_drm=fltarr(m,l)
;
cfrm=fltarr(m,l)
cfthm=fltarr(m,l)
cfphim=fltarr(m,l)
;
 for k=0, l-1 do begin
   dsin_fphi_dthm[*,k]=xder_curl(reform(smooth(sin(theta)*fphim[*,k],2)),theta)       
   d_fr_dthm[*,k]=xder_curl(reform(smooth(frm[*,k],2)),theta)
 endfor
 for j=1,m-1 do begin
   dr_fth_drm[j,*]=xder_curl(smooth(r[*]*fthm[j,*],2),r)   
   dr_fphi_drm[j,*]=xder_curl(smooth(r[*]*fphim[j,*],2),r)
 endfor
 for k=0, l-1 do begin
  for j=1,m-1 do begin
     cfrm[j,k]=rsinth1[j,k]*dsin_fphi_dthm[j,k]
     cfthm[j,k]=rsinth1[j,k]*d_fr_dphim[j,k]-dr_fphi_drm[j,k]/r[k]
     cfphim[j,k]=dr_fth_drm[j,k]/r[k]-d_fr_dthm[j,k]/r[k]
  endfor
 endfor
mhm=cfrm*frm + cfthm*fthm + cfphim*fphim
mhel=(mhtm-mhm)/(4.e-7)

for k=0, l-1 do begin
 mhel[*,k]=mhel[*,k]/rh[k]
endfor
;
tmhel=transpose(mhel)
rgood=0.95
dd=where(r le rgood)
contour,smooth(tmhel[dd,*],2),x[dd,*],y[dd,*],/fi,nl=64,/iso
;
save,filename='mhelicity.sav',r,th,x,y,mhel
;!p.multi=[0,2,1] & !p.charsize=1.3 
;contour,transpose(smooth(khel,3)),x,y,/fi,nl=64,/iso,title='!4a!8=-!4s!8<!4x!8`.u`>/3!6'
;plot,r,smooth(etat,3),/ylog,xtitle='!8r!6',ytitle='!4g!8!Dt!N!6(cm!U!82!N!6s!U!8-1!N!6)!6',$
;       title='!4g!8!Dt!N=!4s!8u!D!6rms!8!U2!N/3!6',thick=3


end
